Volume 2: The Fourier Transform.

<Mark Rose>
<Section 2>
<12/4/18>

Part 1: The Discrete Fourier Transform

In [147]:
from matplotlib import pyplot as plt
from scipy.io import wavfile
import scipy.fftpack
import scipy.signal
import numpy as np
import IPython
import imageio
In [38]:
plt.rcParams["figure.dpi"] = 300             # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3)      # Change plot size / aspect (you may adjust this).
In [83]:
class SoundWave(object):
    """A class for working with digital audio signals."""

    # Problem 1.1
    def __init__(self, rate, samples):
        """Set the SoundWave class attributes.

        Parameters:
            rate (int): The sample rate of the sound.
            samples ((n,) ndarray): NumPy array of samples.
        """
        self.rate = rate                                               #Save the rate as an attribute
        self.samples = samples                                         #Save the samples as an attribute
        return
        raise NotImplementedError("Problem 1.1 Incomplete")

    # Problems 1.1 and 1.7
    def plot(self, force = False):
        """Plot the graph of the sound wave (time versus amplitude)."""
        time = len(self.samples)/self.rate                             #Calculate the time by dividing the length of the samples devided by the rate
        if force == False:                                             #Plot the sound by itself if force is false
            plt.plot(np.linspace(0,time, len(self.samples)), self.samples) #Plot the values according to time and the samples
            plt.xlabel("Time (Seconds)")
            plt.ylabel("Samples")
            plt.title("Sound File Graph")
            plt.ylim((-32768, 32767))                                  #Set the y limits
        else:                                                          #Create two subplots if force is true
            dft = scipy.fftpack.fft(self.samples)                      #Calculate the dft
            freq = np.arange(len(dft))/len(dft)*self.rate              #Calculate the frequencies
            fig, axes= plt.subplots(1,2)
            axes[0].plot(np.linspace(0, time, len(self.samples)), self.samples)  #graph the original sound by doing time by samples
            axes[0].set_xlabel("Time (Seconds)")
            axes[0].set_ylabel("Samples")
            axes[0].set_ylim((-32768,32767))  #Set the y limits
            axes[0].set_title('Sound File Graph')
            axes[1].plot(freq[0:len(freq)//2], np.abs(dft[0:len(freq)//2]))   #Make a second subplot with half of the dft magnitues and frequencies
            axes[1].set_xlabel('Frequency (Hz)')                       #Set the axis labels
            axes[1].set_ylabel('Magnitude')
            axes[1].set_title('Graph of Discrete Fourier Transform')
        return
        raise NotImplementedError("Problem 1.1 Incomplete")

    # Problem 1.2
    def export(self, filename, force=False):
        """Generate a wav file from the sample rate and samples. 
        If the array of samples is not of type np.int16, scale it before exporting.

        Parameters:
            filename (str): The name of the wav file to export the sound to.
        """
        scaled = self.samples
        real = np.real(self.samples)
        if self.samples.dtype != 'int16' or force == True:              #If the samples are not of type int16 or if the boolean force is True, scale the values
            scaled = np.int16((real*32767.0/max(abs(real))))            
        wavfile.write(filename, self.rate, scaled)                      #Write a new sound file with the given filename, same rate, and new scaled values
        return
        raise NotImplementedError("Problem 1.2 Incomplete")
    
    # Problem 1.4
    def __add__(self, other):
        """Combine the samples from two SoundWave objects.

        Parameters:
            other (SoundWave): An object containing the samples to add
                to the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the combined samples.

        Raises:
            ValueError: if the two sample arrays are not the same length.
        """
        if len(self.samples) != len(other.samples):                     #If the lengths are not the same raise an error
            raise ValueError("Cannot add because the sample arrays from the SoundWave objects are not the same size")
        combo = self.samples + other.samples                            #Add the two sounds together
        return SoundWave(self.rate, combo)                              #Return a new SoundWave object with the two sounds and the original rate
        raise NotImplementedError("Problem 1.4 Incomplete")

    # Problem 1.4
    def __rshift__(self, other):
        """Concatentate the samples from two SoundWave objects.

        Parameters:
            other (SoundWave): An object containing the samples to concatenate
                to the samples contained in this object.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        if self.rate != other.rate:                                     #If the rates are not equal, raise an error
            raise ValueError("The two sample rates are not equal")
        combo = np.hstack((self.samples, other.samples))                #Stack the two vectors on top of each other
        return SoundWave(self.rate, combo)                              #Return a sound wave object with the original rate and new combination
        raise NotImplementedError("Problem 1.4 Incomplete")
    
    # Problem 2.1
    def __mul__(self, other):
        """Convolve the samples from two SoundWave objects using circular convolution.
        
        Parameters:
            other (SoundWave): An object containing the samples to convolve
                with the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the convolved samples.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        my_samples = self.samples
        other_samples = other.samples
        if self.rate != other.rate:
            raise ValueError("The two sample rates are not equal")       #If the rates are not equal, raise an error
        if len(self.samples) > len(other.samples):
            num = len(self.samples) - len(other.samples)
            padding = np.zeros(num)
            other_samples = np.hstack((other_samples, padding))
        elif len(other.samples) > self.samples:
            num = len(other.samples) - len(self.samples)
            padding = np.zeros(num)
            my_samples = np.hstack(my_samples, padding)
        conv = scipy.fftpack.ifft(scipy.fftpack.fft(my_samples)*scipy.fftpack.fft(other_samples)).real
        return(SoundWave(self.rate, conv))
            
        raise NotImplementedError("Problem 2.1 Incomplete")

    # Problem 2.2
    def __pow__(self, other):
        """Convolve the samples from two SoundWave objects using linear convolution.
        
        Parameters:
            other (SoundWave): An object containing the samples to convolve
                with the samples contained in this object.
        
        Returns:
            (SoundWave): A new SoundWave instance with the convolved samples.

        Raises:
            ValueError: if the two sample rates are not equal.
        """
        if self.rate != other.rate:
            raise ValueError("The two sample rates are not equal")       #If the rates are not equal, raise an error
        n = len(self.samples)
        m = len(other.samples)
        i=0
        while 2**i < n+m-1:
            i+=1
        pad1 = np.zeros(2**i-n)
        pad2 = np.zeros(2**i-m)
        my_samples = np.hstack((self.samples, pad1))
        other_samples = np.hstack((other.samples, pad2))
        conv = scipy.fftpack.ifft(scipy.fftpack.fft(my_samples)*scipy.fftpack.fft(other_samples)).real
        return(SoundWave(self.rate, conv[0:n+m-1]))
        raise NotImplementedError("Problem 2.2 Incomplete")

    # Problem 2.4
    def clean(self, low_freq, high_freq):
        """Remove a range of frequencies from the samples using the DFT. 

        Parameters:
            low_freq (float): Lower bound of the frequency range to zero out.
            high_freq (float): Higher boound of the frequency range to zero out.
        """
        k_low = int(low_freq*len(self.samples)/self.rate)
        k_high = int(high_freq*len(self.samples)/self.rate)
        n = len(self.samples)
        dft = scipy.fftpack.fft(self.samples)
        dft[k_low:k_high] =0
        dft[n-k_high: n-k_low]=0
        self.samples = scipy.fftpack.ifft(dft).real
        return
        raise NotImplementedError("Problem 2.4 Incomplete")

Problem 1.1

  • Implement SoundWave.__init__().
  • Implement SoundWave.plot().
  • Use the scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.
In [14]:
rate, samples = wavfile.read('tada.wav')                    #Get the rate and samples from the file
tada = SoundWave(rate, samples)                             #Create a SoundWave object
tada.plot()                                                 #Plot the resulting SoundWave

Problem 1.2

  • Implement SoundWave.export().
  • Use the export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).
  • Use IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.
In [15]:
IPython.display.Audio(filename="tada.wav")                #Play the original tada sound
Out[15]:
In [16]:
tada.export("tada_noforce.wav", False)                     #Create a new sound file
IPython.display.Audio(filename="tada_noforce.wav")         #Play the new sound file
Out[16]:
In [17]:
tada.export("tada_withforce.wav", True)                               #Amplify the original tada sound and create a new file
IPython.display.Audio(filename="tada_withforce.wav")                  #Play the amplified sound
Out[17]:

Problem 1.3

  • Implement generate_note().
  • Use generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.
In [18]:
def generate_note(frequency, duration):
    """Generate an instance of the SoundWave class corresponding to 
    the desired soundwave. Uses sample rate of 44100 Hz.
    
    Parameters:
        frequency (float): The frequency of the desired sound.
        duration (float): The length of the desired sound in seconds.
    
    Returns:
        sound (SoundWave): An instance of the SoundWave class.
    """
    sample_space = np.linspace(0,duration, 44100*duration)               #Create a group of values from 0 to the duration
    samples = np.sin(2*np.pi*sample_space*frequency)                     #Use the function to get those values
    rate = 44100                                                         #Use the predefined rate and samples to make a SoundWave object
    note = SoundWave(rate, samples)
    return(note)                                                         #Now return the note
    raise NotImplementedError("Problem 1.3 Incomplete")
In [19]:
my_note = generate_note(440, 2)                                           #Create an A note for 2 seconds
IPython.display.Audio(rate=my_note.rate, data= my_note.samples)           #Play the A note
Out[19]:

Problem 1.4

  • Implement SoundWave.__add__().
  • Generate a three-second A minor chord (A, C, and E) and embed it in the notebook.
  • Implement SoundWave.__rshift__().
  • Generate the arpeggio A$\,\rightarrow\,$C$\,\rightarrow\,$E, where each tone lasts one second, and embed it in the notebook.
In [20]:
a_note = generate_note(440,3)                                                   #Create an A, C, and E sounds that all last 3 seconds 
c_note = generate_note(523.25,3)
e_note = generate_note(659.25,3)
a_minor_chord = a_note+c_note+e_note                                            #Create an A minor chord
IPython.display.Audio(rate=a_minor_chord.rate, data = a_minor_chord.samples)    #Play the chord
Out[20]:
In [21]:
a_note = generate_note(440,1)                                                        #Make an a, c, and e note that are all one second long
c_note = generate_note(523.25,1)
e_note = generate_note(659.25,1)
a_minor_arpeggio = a_note >> c_note >> e_note                                        #Concatenate the three sounds
IPython.display.Audio(rate=a_minor_arpeggio.rate, data = a_minor_arpeggio.samples)   #Play the three sounds one after another
Out[21]:

Problem 1.5

  • Implement simple_dft() with the formula for $c_k$ given below.
  • Use np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).
$$ c_k = \frac{1}{n}\sum_{\ell=0}^{n-1} f_\ell e^{-2 \pi i k \ell\, /\, n} $$
In [22]:
def simple_dft(samples):
    """Compute the DFT of an array of samples.

    Parameters:
        samples ((n,) ndarray): an array of samples.
    
    Returns:
        ((n,) ndarray): The DFT of the given array.
    """
    n = len(samples)
    matrix_vals = np.reshape(np.arange(n), (n,1)) * np.reshape(np.arange(n), (1,n))  #Broadcast to make a matrix for the l and k values
    dft = (samples * np.exp(-2*np.pi*1j*matrix_vals/n)).sum(axis=1)/n                #Use the above equation to calculate the dft
    return(dft)
    raise NotImplementedError("Problem 1.5 Incomplete")
In [23]:
for n in range(6,20):
    samples = np.random.random(n)                    #Make a bunch of random samples
    my_dft = simple_dft(samples)                     #Check my dft vs their dft
    given_dft = scipy.fftpack.fft(samples)          
    print(np.allclose(my_dft*n, given_dft))          #Print the results

True
True
True
True
True
True
True
True
True
True
True
True
True
True

Problem 1.6

  • Implement simple_fft().
  • Generate an array of $8192$ random samples and take its DFT using simple_dft(), simple_fft(), and scipy.fftpack.fft(). Print the runtimes of each computation.
  • Use np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).
In [24]:
def simple_fft(samples, threshold=1):
    """Compute the DFT using the FFT algorithm.
    
    Parameters:
        samples ((n,) ndarray): an array of samples.
        threshold (int): when a subarray of samples has fewer
            elements than this integer, use simple_dft() to
            compute the DFT of that subarray.
    
    Returns:
        ((n,) ndarray): The DFT of the given array.
    """
    def split(g,N):
        n = len(g)
        if n <= N:
            return(n*simple_dft(g))                                    #Use the predifined function with a small g
        else:
            even = split(g[::2],N)
            odd = split(g[1::2],N)                                     #Split into evens and odds
            z = np.zeros(n)
            z = np.exp(-2*np.pi*1j*np.arange(n)/n)                    #calculate the exponential parts using broadcasting
            m = n // 2
            return(np.hstack((even + z[:m] * odd, even + z[m:]*odd))) #return the concatenation of the two arrays
    return(split(samples, threshold)/len(samples))                    #Start the recursion
    raise NotImplementedError("Problem 1.6 Incomplete")
In [25]:
powers_of_two = [1,2,4,8,16,32,64,128]                  #Create a list of powers of two
for n in powers_of_two:
    samples = np.random.random(n)                       #Generate random sample vectors
    my_dft = simple_fft(samples)
    given_dft = scipy.fftpack.fft(samples)              #Calculate the samples using my function, and the predifined one
    print(np.allclose(my_dft*n, given_dft))             #Check to see if the functions have the same result

True
True
True
True
True
True
True
True
In [26]:
samples = np.random.random(8192)                             
%time simple_dft(samples)                                 #This line times the simple dft algorithm                
%time simple_fft(samples)                                 #This line times the simple fft algorithm
%time scipy.fftpack.fft(samples)                          #This line times the predifined dft algorithm

CPU times: user 4.01 s, sys: 733 ms, total: 4.74 s
Wall time: 4.74 s
CPU times: user 197 ms, sys: 26 µs, total: 197 ms
Wall time: 196 ms
CPU times: user 210 µs, sys: 23 µs, total: 233 µs
Wall time: 221 µs
Out[26]:
array([ 4.13238298e+03 +0.j        ,  5.52243104e+00-10.20770873j,
       -2.55857247e+00+34.23527382j, ..., -3.89716644e+01-18.7670011j ,
       -2.55857247e+00-34.23527382j,  5.52243104e+00+10.20770873j])

Problem 1.7

  • Modify SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.
  • Display the plot of the DFT of an A tone.
  • Display the plot of the DFT of an A minor chord.
In [27]:
a_note = generate_note(440,2)                                       #Create an A note
a_note.plot(True)                                                   #Plot that note
In [28]:
a_note = generate_note(440,2)                                              #The graph on here and the one in the book 
c_note = generate_note(523.25,2)                                           #match when time = 2 seconds
e_note = generate_note(659.25,2)
a_minor_chord = a_note+c_note+e_note                                       #Create a chord with the notes a, c, and e
a_minor_chord.plot(True)                                                   #Use the plotting method to plot the sounds

Problem 1.8

Use the DFT to determine the individual notes that are present in mystery_chord.wav.

In [29]:
mrate, msamples = wavfile.read('mystery_chord.wav')                            #Get the rate and the samples from the sound file
myst= SoundWave(mrate, msamples)                                               #Create a SoundWave object
dft=scipy.fftpack.fft(msamples)
dft=np.abs(dft)                                                                #Get the dft and the magnitudes
freq = np.arange(len(dft))/len(dft)*mrate                                      #calculate the frequencies
vals = np.argsort(dft[0:len(freq)//2])[-4:]                                    #Take the largest 3 from the first half
sounds = freq[vals]
In [30]:
sounds                                                                         #These are the frequencies of the sounds
Out[30]:
array([587.5 , 523.25, 784.  , 440.  ])

The notes are D, C, G, and A

Part 2: Convolution and Filtering.

Problem 2.1

  • Implement SoundWave.__mul__() for circular convolution.
  • Generate 2 seconds of white noise at the same sample rate as tada.wav.
  • Compute the circular convolution of tada.wav and the white noise. Embed the result in the notebook.
  • Append the circular convolution to itself and embed the result in the notebook.
In [34]:
rate, samples = wavfile.read('tada.wav')                            
tada = SoundWave(rate, samples)                                     #create a SoundWave object with the rate and samples of tada
white_noise = SoundWave(rate, np.random.randint(-32767,32767, rate*2, dtype=np.int16))  #Create some random white noise in a SoundWave object
my_noise = white_noise * tada                                       #Convolute the two noises
print(my_noise.rate)
my_noise = my_noise >> my_noise                                     
IPython.display.Audio(rate=my_noise.rate, data = my_noise.samples)  #Give the sound

22050
Out[34]:

Problem 2.2

  • Implement SoundWave.__pow__() for linear convolution.
  • Time the linear convolution of CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().
  • Embed the two original sounds and their convolutions in the notebook. Check that the convolutions with SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.
In [65]:
cgc_rate, cgc_samples = wavfile.read('CGC.wav')
cgc = SoundWave(cgc_rate, cgc_samples)
gcg_rate, gcg_samples = wavfile.read('GCG.wav')
gcg = SoundWave(gcg_rate, gcg_samples)

%time new_sound = cgc**gcg                                                #This line times the simple dft algorithm                
%time defined_samples = scipy.signal.fftconvolve(cgc_samples, gcg_samples)  #This line times the simple fft algorithm

CPU times: user 475 ms, sys: 12.3 ms, total: 488 ms
Wall time: 82.2 ms
CPU times: user 65.7 ms, sys: 167 µs, total: 65.9 ms
Wall time: 11.6 ms
In [66]:
IPython.display.Audio(rate=new_sound.rate, data = new_sound.samples)
Out[66]:
In [67]:
IPython.display.Audio(rate=cgc_rate, data = defined_samples)
Out[67]:

Problem 2.3

Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav. Embed the two original sounds and their convolution in the notebook.

In [92]:
c_rate, c_samples = wavfile.read('chopin.wav')
chop = SoundWave(c_rate, c_samples)
b_rate, b_samples = wavfile.read('balloon.wav')
ball = SoundWave(b_rate, b_samples)
chopin_echo = chop ** ball
In [69]:
IPython.display.Audio(rate=c_rate, data = c_samples)
Out[69]:
In [70]:
IPython.display.Audio(rate=b_rate, data = b_samples)
Out[70]:
In [71]:
IPython.display.Audio(rate=chopin_echo.rate, data = chopin_echo.samples)
Out[71]:

Problem 2.4

  • Implement SoundWave.clean().
  • Clean noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.
  • Clean noisy2.wav. Embed the original and the cleaned versions in the notebook.
In [122]:
n1_rate, n1_samples = wavfile.read('noisy1.wav')
noisy1 = SoundWave(n1_rate, n1_samples)
n2_rate, n2_samples = wavfile.read('noisy2.wav')
noisy2 = SoundWave(n2_rate, n2_samples)
In [123]:
IPython.display.Audio(rate=noisy1.rate, data = noisy1.samples)
Out[123]:
In [124]:
noisy1.clean(1250,2600)
IPython.display.Audio(rate=noisy1.rate, data = noisy1.samples)
Out[124]:
In [125]:
IPython.display.Audio(rate=noisy2.rate, data = noisy2.samples)
Out[125]:
In [126]:
noisy2.plot(force=True)
In [127]:
noisy2.clean(1300,4400)                                         #I determined from the above graph to filter from 1300 to 4400
IPython.display.Audio(rate=noisy2.rate, data = noisy2.samples)
Out[127]:

Problem 2.5

  • Clean vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.
  • Recombine the left and right channels and embed the result in the notebook.
In [315]:
vuvu_rate, vuvu_samples = wavfile.read('vuvuzela.wav')
vuvu_left = SoundWave(vuvu_rate, vuvu_samples[:,0])
vuvu_right = SoundWave(vuvu_rate, vuvu_samples[:,1])
vuvu_left.clean(200,500)
vuvu_right.clean(200,500)
new_sample = np.vstack((vuvu_left.samples, vuvu_right.samples))
print("Original Signal")

Original Signal
In [316]:
IPython.display.Audio(rate=vuvu_rate, data =vuvu_samples.T )
Out[316]:
In [317]:
print('New Signal')

New Signal
In [318]:
IPython.display.Audio(rate=vuvu_rate, data =new_sample )
Out[318]:

Problem 2.6

  • Clean up license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.
  • Display the original and cleaned images.
In [337]:
im = imageio.imread('license_plate.png')
fft2 = scipy.fftpack.fft2(im)
plt.imshow(np.log(np.abs(fft2)),cmap='gray')
plt.show()
In [338]:
mask = np.log(np.abs(fft2)) >12
fft2[mask] = np.mean(fft2)
plt.imshow(np.log(np.abs(fft2)),cmap='gray')
plt.show()
In [339]:
new_im = scipy.fftpack.ifft2(fft2).real
plt.imshow(new_im, cmap='gray')
Out[339]:
<matplotlib.image.AxesImage at 0x7ff7d40bf1d0>

The year on the sticker is 2013